       PROGRAM  LSSPHRD                                                 00010000
CCC                                                                     00020000
CCCCC -------LIGHT SCATTERING BY AN ARBITRALLY ORIENTED SPHEROID--------00030000
CCC          INTENSITY MATRIX FOR AN INCIDENCE OF NATURAL OR            00040000
CCC          LINEARLY POLARIZED LIGHT.                                  00050000
CCC             THE STATE OF POLARIZATION FOR THE INCIDENT LIGHT        00060000
CCC             IS SPECIFIED IN THE PHI = 0 PLANE.                      00070000
CCC             THE SCATTERED FIELD IS BIVEN IN THE ABSOLUTE COORDINATE 00080000
CCC             SYSTEM RELATED TO THE DIRECTION OF INCIDENT LIGHT.      00090000
CCC             NSHAPE = 1 :  PROLATE SPHEROID(IPO= 1)                  00100000
CCC             NSHAPE = 2 :  OBLATE  SPHEROID(IPO=-1)                  00110000
CCC             ZETA = INCIDENCE ANGLEG                                 00120000
CCC             RFR = COMPLEX REFRACTIVE INDEX (RFRRL, RFRIM)           00130000
CCC             A = SEMI-MAJOR AXIS ; B = SEMI-MINOR AXIS ; ABYB = A/B  00140000
CCC             NMAX .LE. NMX ;   MMAX .LE. MMX ;   INCANG.LE. NINC     00150034
CCC             SIMPSON LAW IS USED FOR INTEGRATIONS                    00160000
CCC          ORIGINALLY PROGRAMMED IN 1972                              00170000
CCC          ICALC < 0 : SKIP INTENSITY CALCULATION
CCC          REFERENCES: S. ASANO & G. YAMAMOTO (1975), S. ASANO (1979)
CCCCC ----------------- REVISED IN FORTRAN77 IN 1978 ------------------ 00180000
      IMPLICIT  REAL*8(A-H,O-Z)                                         00190000
      PARAMETER  (JMX=150, NMX= 70, MMX=31, NSYS=45, NINC=10)           00191038
      COMPLEX*16 C2,CI,CJ,RFR,RI2,DRI2,CLAMD,C2XI,DUMY1,DUMY2,          00200000
     1           DC(JMX),DI(NMX,JMX),SPCI(NMX),RE3(NMX),DRE3(NMX),      00210036
     2           RI1(NMX),DRI1(NMX),UE1(NMX,NMX),UE3(NMX,NMX),          00220034
     3           UI1(NMX,NMX),VE1(NMX,NMX),VE3(NMX,NMX),VI1(NMX,NMX),   00221034
     4           XE1(NMX,NMX),XE3(NMX,NMX),XI1(NMX,NMX),YE1(NMX,NMX),   00230034
     5           YE3(NMX,NMX),YI1(NMX,NMX),F(NMX,NSYS),G(NMX,NSYS),     00240034
     6           ALP(NMX),BET(NMX),ALP1(NMX,MMX,2),BET1(NMX,MMX,2),     00250034
     7           DATA(2*NMX),TE1,TE2,TM1,TM2,A1,A2,A3,A4                00260034
      COMPLEX*16 CIN(4)/(1.D0,0.D0),(0.D0,-1.D0),(-1.D0,0.D0),          00270000
     +                  (0.D0,1.D0)/                                    00280034
      COMPLEX*16 CJN(4)/(1.D0,0.D0),(0.D0,1.D0),(-1.D0,0.D0),           00290000
     +                  (0.D0,-1.D0)/                                   00300034
      DIMENSION PL(JMX),DPL(JMX),DR(JMX),DE(NMX,JMX),SIGME(NMX),        00310036
     1          SUMD(NMX),SPCE(NMX),RE1(NMX),RE2(NMX),DRE1(NMX),        00320034
     2          JWE(NMX),JWI(NMX),NUMB(MMX),CONV(MMX),SIG0(NMX,NINC),   00330034
     3          CHI0(NMX,NINC),JWM(MMX,NMX),DEM(MMX,NMX,65),PFACT(MMX), 00331034
     4          BSCA(NMX,NMX),EXT(NINC,2),SCA(NINC,2),NAZM(5),          00340034
     5          SANG(91),RLTHT(91),WTHT(91),RLPHI(37),WPHI(37),         00350034
     5          DZETA(10),ZETA(10),PG(10),R(91,37,2,2),P(91,2,2)        00360000
      COMMON  /XYUV/SPCI,RI1,DRI1,RE3,DRE3,DE,RE1,DRE1,SPCE             00370034
      COMMON  /CMPLX/DI                                                 00371034
      COMMON  /BNDRY/UE1,UE3,UI1,VE1,VE3,VI1,XE1,XE3,XI1,YE1,YE3,YI1,   00380000
     +               F,G,RFR                                            00390034
C*                                                                      00400000
      DATA PI,CI,CJ,EPS/3.141592653589793D0,(0.D0,1.D0),(1.D0,0.D0),    00410000
     +                  1.D-10/                                         00420000
      LOGICAL*1 SHAPE                                                   00430000
C*      CALL  ERRSET(207,256,-1,0,0,208)                                00450000
C                                                                       00460000
         DEFINE FILE  9(1000,400,U,IV1)                                 00461034
         NUNIT=9                                                        00462034
      ICALC= 1                                                          00470031
C*    ICALC= -1     ( SKIP OF INTENSITY CALC.)   
      MINT=0                                                            00500000
      MMAX= MMX-1                                                       00510038
      NLMT= NMX                                                         00520038
      STRD=PI/180.D0                                                    00530000
      OMG=1.D0-EPS                                                      00540000
C*      CALL DATE(QDATE)                                                00550034
CC        INPUT DATA OF INCREMENTS OF ANGLES
      READ(5,500) DINC,DSANG,DPHI                                       00560000
  500 FORMAT(' DINC=',F5.0,' DSANG=',F5.0,' DPHI=',F5.0)                00570000
C         SETTING OF THE INCIDENCE ANGLES                               00580000
      INCD=(90.0+EPS)/DINC                                              00590000
      INCANG=INCD+1                                                     00600000
      DO 10 I=1,INCANG                                                  00610000
      DZETA(I)=DINC*(I-1)                                               00620000
   10 ZETA(I)=STRD*DZETA(I)                                             00630000
C         OBSERVATION POINTS (RLTHT,RLPHI)                              00640000
C          SIMPSON LAW FOR INTEGRATION OVER AZIMUTH ANGLES              00650000
      NTHT=(180.0+EPS)/DSANG+1                                          00660000
      SIMPSN=STRD*DSANG/3.D0                                            00670000
      WTHT(1)=SIMPSN                                                    00680000
      WTHT(NTHT)=SIMPSN                                                 00690000
      DO 20 J=1,NTHT                                                    00700000
      SANG(J)=DSANG*(J-1)                                               00710000
      RLTHT(J)=STRD*SANG(J)                                             00720000
      IF(J.EQ.1.OR.J.EQ.NTHT) GO TO 20                                  00730000
      WTHT(J)=2*SIMPSN                                                  00740000
      IF(MOD(J,2).EQ.0) WTHT(J)=4*SIMPSN                                00750000
   20 CONTINUE                                                          00760000
      NPHI=(180.0+EPS)/DPHI                                             00770000
      NPHI=2*(NPHI/2)                                                   00780000
      DPHI=180.0/DFLOAT(NPHI)                                           00790000
      SIMPSN=STRD*DPHI/3.D0/PI                                          00800000
      NPHI1=NPHI+1                                                      00810000
      WPHI(1)=SIMPSN                                                    00820000
      WPHI(NPHI1)=SIMPSN                                                00830000
      DO 30 K=1,NPHI1                                                   00840000
      RLPHI(K)=STRD*DPHI*(K-1)                                          00850000
      IF(K.EQ.1.OR.K.EQ.NPHI1) GO TO 30                                 00860000
      WPHI(K)=2*SIMPSN                                                  00870000
      IF(MOD(K,2).EQ.0) WPHI(K)=4*SIMPSN                                00880000
   30 CONTINUE                                                          00890000
      NAZM(1)=1                                                         00900000
      NAZM(2)=NPHI/4+1                                                  00910000
      NAZM(3)=NPHI/2+1                                                  00920000
      NAZM(4)=3*(NPHI/4)+1                                              00930000
      NAZM(5)=NPHI1                                                     00940000
CC         SPHEROID  **  SIZE P.=C, A/B=ABYB, REFRACT.IND.=RFR          00950000
  400 READ(5,200,END=300) NSHAPE,ABYB,A,RFRRL,RFRIM                     00960000
  200 FORMAT(' NSHAPE=',I2,' ABYB=',F5.0,' A=',F5.0,' RFRRL=',F6.0,     00970000
     +       ' RFRIM=',F7.0 )                                           00980000
      RFR=DCMPLX(RFRRL,RFRIM)                                           00990000
C      WRITE(6,666) QDATE                                               01000000
C 666  FORMAT(1H1,100X,'DATE =', A8)                                    01010000
C      CALL CLOCK                                                       01020000
CC               PROLATE SPHEROID                                       01021034
      SHAPE=NSHAPE.EQ.1                                                 01030034
        IF(SHAPE)  THEN                                                 01040034
      IPO= 1                                                            01060000
      XI=ABYB/DSQRT(ABYB*ABYB-1.D0)                                     01070000
      XI2=XI*XI-1.D0                                                    01080000
      ECCEN=1.D0/XI                                                     01090000
      C=A*ECCEN                                                         01100000
      B=C*DSQRT(XI2)                                                    01110000
      RV=(A*B*B)**(1.D0/3.D0)                                           01120000
      RTAB=DSQRT(A*A-B*B)                                               01130000
      SAREA=0.5D0*PI*(B*B+A*A*B*DARCOS(1.D0/ABYB)/RTAB)                 01140000
      RS=DSQRT(SAREA/PI)                                                01150000
      WRITE(6,111)                                                      01160000
  111 FORMAT(1H0/20X,'*****   LIGHT SCATTERING BY AN ARBITRALLY ORIENTED01170000
     + PROLATE SPHEROIDAL PARTICLE   *****' )                           01180000
        ELSE                                                            01190034
CC              OBLATE SPHEROID                                         01200000
      IPO=-1                                                            01210034
      XI=1.D0/DSQRT(ABYB*ABYB-1.D0)                                     01220000
      XI2=XI*XI+1.D0                                                    01230000
      ECCEN=1.D0/DSQRT(XI2)                                             01240000
      C=A*ECCEN                                                         01250000
      B=C*XI                                                            01260000
      RV=(A*A*B)**(1.D0/3.D0)                                           01270000
      RTAB=DSQRT(A*A-B*B)                                               01280000
      SAREA=0.5D0*PI*(A*A+0.5*A*B*B/RTAB*DLOG((A+RTAB)/(A-RTAB)))       01290000
      RS=DSQRT(SAREA/PI)                                                01300000
      WRITE(6,222)                                                      01310000
  222 FORMAT(1H0/20X,'*****   LIGHT SCATTERING BY AN ARBITRALLY ORIENTED01320000
     + OBLATE SPHEROIDAL PARTICLE   *****' )                            01330000
        END IF                                                          01340034
      C1=C                                                              01341034
      C2=C1*RFR                                                         01350000
      C3=C2                                                             01360000
C
      WRITE(6,100) C,RFR,ABYB,ECCEN,A,RV,B,RS                           01370000
  100 FORMAT(1H0//28X,'C   =',F7.4,17X,'REFRACTIVE INDEX =('            01380000
     + ,F6.3,',',F7.4,')'/28X,'A/B =',F7.3,17X,'ECCENTRICITY =',F8.5/   01390000
     + 28X,'K*A =',F7.3,17X,'EQUI-RADIUS  K*RV =',F8.5/28X,'K*B =',F7.3,01400000
     + 17X,'EQUI-RADIUS  K*RG =',F8.5 )                                 01410000
      WRITE(6,106) DINC,DSANG,DPHI                                      01420000
  106 FORMAT(1H0/28X,'DINCANG =',F6.2,'DEG',12X,'DSANG =',F6.2,'DEG',   01430000
     + 12X,'DPHI =',F6.2,'DEG'//)                                       01440000
C
CC           EMPIRICAL NUMBER OF EXPANSION TERMS
      IF(ABYB.GE.3.0) GO TO 40                                          01450000
      KMIN0=0.8*A+20                                                    01460000
      IF(A.LE.1.0) KMIN0=14                                             01470000
      GO TO 50                                                          01480000
   40 KMIN0=0.8*A+30                                                    01490000
      IF(A.LT.5.0) KMIN0=0.8*A+24                                       01500000
   50 CONTINUE                                                          01510000
      MAXJ1=40+KMIN0/2+3.0*C1                                           01520000
      MAXJ2=40+KMIN0/2+3.0*CDABS(C2)                                    01530000
       IF(MAXJ1.GE.JMX)  MAXJ1= JMX-2                                   01540036
       IF(MAXJ2.GE.JMX)  MAXJ2= JMX-2                                   01550036
CC           INCIDENCE ANGLES -ZETA- FROM (ELV,BETA)                    01560000
      BIG=0.D0                                                          01570000
      DO 9000 I=1,INCANG                                                01580000
      CINC= DCOS(ZETA(I))                                               01590034
      SINC= DSIN(ZETA(I))                                               01600034
      C2INC=CINC*CINC                                                   01600134
      S2INC=SINC*SINC                                                   01601034
        IF(SHAPE)  THEN                                                 01610034
      PG(I)=PI*B*B*DSQRT(ABYB*ABYB*S2INC+C2INC)                         01620000
        ELSE                                                            01630034
      PG(I)=PI*A*A*DSQRT(S2INC/ABYB/ABYB+C2INC)                         01640034
        END IF                                                          01641034
      DO 9010 L=1,2                                                     01650034
      EXT(I,L)=0.                                                       01660000
 9010 SCA(I,L)=0.                                                       01670000
      CONV(I)=0.D0                                                      01680000
      IF(BIG-ZETA(I)) 9020,9000,9000                                    01690000
 9020 IZETA=I                                                           01700000
      BIG=ZETA(IZETA)                                                   01710000
 9000 CONTINUE                                                          01720000
CC          AZIMUTHAL DEPENDENCE  ***  M                                01730000
      INDEX=0                                                           01740000
      M=0                                                               01750000
 9050 CRIT0=EXT(IZETA,1)+EXT(IZETA,2)                                   01760000
      MM=M+1                                                            01770000
      M2=M+M                                                            01780000
      WRITE(6,102) M                                                    01790000
  102 FORMAT(1H /11X,'*   M =',I2 )                                     01800000
      KMIN=KMIN0-DFLOAT(M)                                              01810000
C           CORRECTION FOR THE (M=0) TERM CONVERGENCE                   01820000
      IF(M.EQ.0.AND.A.GT.14. ) KMIN=KMIN0+10                            01830000
      IF(KMIN.LE.10) KMIN=10                                            01840000
      IF(KMIN.GE.NLMT) KMIN=NLMT                                        01850000
      KMAX=KMIN                                                         01860000
      PFACT(MM)=1.D0                                                    01870000
      MINDX=0                                                           01880033
    8 MINDX= MINDX+1                                                    01881033
      PFACT(MM)=PFACT(MM)*(2*MINDX-1)                                   01890033
      IF(MINDX.LT.M)  GO TO 8                                           01891033
CC        SPHEROIDAL RADIAL FUNCTIONS   **   SUBROUTINE                 01900000
      DO 9060 I=1,JMX                                                   01920036
      DO 9060 N=1,KMAX                                                  01921034
      DE(N,I)= 0.D0                                                     01930000
 9060 DI(N,I)= (0.D0,0.D0)                                              01940034
      C1XI=-(C1*XI)**2+IPO*M*M/XI2                                      01950000
      DO 1060 NN=1,KMAX                                                 01960000
      N=M+NN-1                                                          01970000
      L=N-M                                                             01980000
       IF(SHAPE) THEN                                                   01990000
        CALL  RPRSWF(C1,XI,M,N,MINT,1,KMAX,MAXJ1,JW,DR,ALAMD,SIGME(NN), 02000031
     1               SUMD(NN),RE1(NN),DRE1(NN),RE2(NN),DRE2)            02010031
       ELSE                                                             02020034
        CALL  ROBSWF(C1,XI,M,N,MINT,1,KMAX,MAXJ1,JW,DR,ALAMD,SIGME(NN), 02030034
     1               SUMD(NN),RE1(NN),DRE1(NN),RE2(NN),DRE2)            02040031
       END IF                                                           02050034
      SPCE(NN)=ALAMD+C1XI                                               02051034
      JWE(NN)=JW                                                        02060000
      IF(MOD(L,2).EQ.0) JWE(NN)=JW+1                                    02070000
      JJA=JWE(NN)                                                       02080000
      JWM(MM,NN)=JJA                                                    02090000
      DO 1065 J=1,JJA,2                                                 02100000
      JP=(J/2)+1                                                        02110000
      DEM(MM,NN,JP)=DR(J)                                               02120000
 1065 DE(NN,J)=DR(J)                                                    02130000
      RE3(NN)=RE1(NN)+CI*RE2(NN)                                        02140000
      DRE3(NN)=DRE1(NN)+CI*DRE2                                         02150000
      NK=NN                                                             02160000
      IF(DABS(RE1(NN)).LT.1.D-70.OR.DABS(DRE2).GT.1.D+70) GO TO 1040    02170000
 1060 CONTINUE                                                          02180000
      GO TO 1050                                                        02190000
 1040 KMAX=NK                                                           02200000
      IF(KMIN.GE.KMAX) KMIN=KMAX                                        02210000
 1050 CONTINUE                                                          02220000
      C2XI=-(C2*XI)**2+IPO*M*M/XI2                                      02230000
      C3XI=-(C3*XI)**2+IPO*M*M/XI2                                      02240000
      DO 1070 NN=1,KMAX                                                 02250000
      N=M+NN-1                                                          02260000
      L=N-M                                                             02270000
      IF(RFRIM.LT.1.D-9) GO TO 1030                                     02280000
       IF(SHAPE)  THEN                                                  02290034
        CALL  CPRSWF(C2,XI,M,N,MINT,2,KMAX,MAXJ2,JW,DC,CLAMD,DUMY1,     02300032
     +               DUMY2,RI1(NN),DRI1(NN),RI2,DRI2)                   02310032
       ELSE                                                             02320034
        CALL  COBSWF(C2,XI,M,N,MINT,2,KMAX,MAXJ2,JW,DC,CLAMD,DUMY1,     02330034
     +              DUMY2,RI1(NN),DRI1(NN),RI2,DRI2)                    02340032
       END IF                                                           02350034
      SPCI(NN)=CLAMD+C2XI                                               02351034
      JWI(NN)=JW                                                        02360000
      IF(MOD(L,2).EQ.0) JWI(NN)=JW+1                                    02370000
      JJB=JWI(NN)                                                       02380000
      DO 1075 J=1,JJB,2                                                 02390000
 1075 DI(NN,J)=DC(J)                                                    02400000
      GO TO 1070                                                        02410000
 1030  IF(SHAPE)  THEN                                                  02420034
        CALL  RPRSWF(C3,XI,M,N,MINT,2,KMAX,MAXJ2,JW,DR,ALAMD,DUMY3,     02430032
     +               DUMY4,R1,DR1,R2,DR2)                               02440032
       ELSE                                                             02450034
        CALL  ROBSWF(C3,XI,M,N,MINT,2,KMAX,MAXJ2,JW,DR,ALAMD,DUMY3,     02460034
     +               DUMY4,R1,DR1,R2,DR2)                               02470032
       END IF                                                           02471034
      SPCI(NN)=CJ*(ALAMD+C3XI)                                          02480034
      RI1(NN)=CJ*R1                                                     02490000
      DRI1(NN)=CJ*DR1                                                   02500000
      JWI(NN)=JW                                                        02510000
      IF(MOD(L,2).EQ.0) JWI(NN)=JW+1                                    02520000
      JJB=JWI(NN)                                                       02530000
      DO 1035 J=1,JJB,2                                                 02540000
 1035 DI(NN,J)=CJ*DR(J)                                                 02550000
 1070 CONTINUE                                                          02560000
      DO 1080 NN=1,KMAX                                                 02570000
      JEN=JWE(NN)                                                       02580000
      JIN=JWI(NN)                                                       02590000
      DO 1090 LL=1,KMAX                                                 02600000
CC        INTEGRALS A-I  &  PARAMETERS U-Y                              02610000
        CALL  UVXYAI(C1,C2,XI,M,JEN,JIN,NN,LL,UE1(LL,NN),UE3(LL,NN),    02620000
     +      UI1(LL,NN),VE1(LL,NN),VE3(LL,NN),VI1(LL,NN),XE1(LL,NN),     02630034
     +      XE3(LL,NN),XI1(LL,NN),YE1(LL,NN),YE3(LL,NN),YI1(LL,NN),IPO) 02640034
CC        PRECALCULATIONS FOR SCATTERING CROSS SECTION                  02650000
      BSCA(NN,LL)=0.D0                                                  02660000
      NAB=IABS(NN-LL)                                                   02670000
      IF(MOD(NAB,2).EQ.1) GO TO 1090                                    02680000
      DO 7020 IA=1,JEN,2                                                02690000
      I=IA-MOD(NN,2)                                                    02700000
      Q= DE(NN,IA)                                                      02710000
      IF(M.EQ.0) GO TO 7050                                             02720000
      DO 7060 IC=1,M2                                                   02730000
 7060 Q=Q*DFLOAT(I+IC)                                                  02740000
 7050 Q=2.0*Q*DFLOAT((M+I)*(M+I+1))/DFLOAT(M2+I+I+1)                    02750000
      IF(M.EQ.0) Q=2.0*Q                                                02760000
      DEAB=Q*DE(LL,IA)                                                  02770000
      IF(I.GT.NN.AND.DABS(DEAB).LT.1.D-30) GO TO 1090                   02780000
      BSCA(NN,LL)=BSCA(NN,LL)+DEAB                                      02790000
 7020 CONTINUE                                                          02800000
 1090 CONTINUE                                                          02810000
 1080 CONTINUE                                                          02820000
CC            OBLIQUELY INCIDENT WAVES                                  02830000
      DO 3010 IN=1,INCANG                                               02840000
      ZETA0=ZETA(IN)                                                    02850000
      IF(ZETA0.LE.1.D-8) ZETA0=1.D-8                                    02860000
      ETA0=DCOS(ZETA0)                                                  02870000
      ETA20=1.0-ETA0*ETA0                                               02880000
      SQRE0=DSQRT(ETA20)                                                02890000
      PL(1)=PFACT(MM)*(SQRE0**M)                                        02900000
      PL(2)=(2*M+1)*ETA0*PL(1)                                          02910000
      DO 4030 I=1,MAXJ1                                                 02920000
      MI=M+I                                                            02930000
      FI1=I+1                                                           02940000
      PL(I+2)=((2*MI+1)*ETA0*PL(I+1)-(MI+M)*PL(I))/FI1                  02950000
 4030 DPL(I)=(MI*ETA0*PL(I)-(MI-M)*PL(I+1))/ETA20                       02960000
      DO 3040 NN=1,KMAX                                                 02970000
      N=M+NN-1                                                          02980000
      NM=N-M                                                            02990000
      JEN=JWE(NN)                                                       03000000
      NMODN=MOD(N,4)+1                                                  03010000
      FA=0.D0                                                           03020000
      GA=0.D0                                                           03030000
      S0=0.D0                                                           03040000
      DS0=0.D0                                                          03050000
      DO 3000 IR=1,JEN,2                                                03060000
      I=IR+MOD(NM,2)-1                                                  03070000
      DEPL0=DE(NN,IR)*PL(I+1)                                           03080000
      S0=S0+DEPL0                                                       03090000
      DS0=DS0+DE(NN,IR)*DPL(I+1)                                        03100000
      IF(M+I.EQ.0) GO TO 3020                                           03110000
      Q1=DE(NN,IR)/DFLOAT((M+I)*(M+I+1))                                03120000
      FA=FA+Q1*PL(I+1)                                                  03130000
      GA=GA-SQRE0*Q1*DPL(I+1)                                           03140000
      IF(I.GE.20.AND.DABS(DEPL0).LT.1.D-20) GO TO 3080                  03150000
      GO TO 3000                                                        03160000
 3020 GA= DE(NN,1)*(1.D0-ETA0)/SQRE0                                    03170000
      FA=0.D0                                                           03180000
 3000 CONTINUE                                                          03190000
 3080 Q2=4.D0                                                           03200000
      IF(M.EQ.0) Q2=2.D0                                                03210000
      F(NN,IN)=4.D0*M*CJN(NMODN)*FA/(SIGME(NN)*SQRE0)                   03220000
      G(NN,IN)=Q2*CJN(NMODN)*GA/SIGME(NN)                               03230000
      SIG0(NN,IN)=DFLOAT(M)*S0/SQRE0                                    03240000
      CHI0(NN,IN)=-SQRE0*DS0                                            03250000
 3040 CONTINUE                                                          03260000
 3010 CONTINUE                                                          03270000
CC           DETERMINATION OF EXPANSION COEFFICIENTS  -  ALP & BET      03280000
      NMAX=2*(KMIN/2)                                                   03290000
      NUMB(MM)=NMAX                                                     03300000
      DO 5000 IMODE=1,2                                                 03310000
      IF(IMODE.EQ.1) WRITE(6,103)                                       03320000
  103 FORMAT(1H ,14X,'A)  TM-MODE POLARIZATION ' )                      03330000
C             THE BOUNDARY CONDITIONS                                   03340000
        CALL BOUNDR(INCANG,M,KMIN,IMODE,RE1,RE2,RI1)                    03350000
      DO 3050 IN=1,INCANG                                               03360000
      DO 3030 NE=1,NMAX                                                 03370000
      N=M+NE-1                                                          03380000
      NMODN=MOD(N,4)+1                                                  03390000
      ALP(NE)= CIN(NMODN)*DI(IN,NE)                                     03400000
      BET(NE)= CIN(NMODN)*DI(IN,NE+NMX)                                 03410034
      DATA(NE)=ALP(NE)                                                  03420000
      DATA(NE+NMX)=BET(NE)                                              03430034
 3030 CONTINUE                                                          03440000
       INDEX=INDEX+1                                                    03450000
       WRITE(NUNIT'INDEX) DATA                                          03460000
CC            SCATTERING CROSS-SECTIONS                                 03470000
      IF(ZETA(IN).LT.1.D-8.AND.M.NE.1) GO TO 3050                       03480000
      IF(M.GE.3.AND.CONV(IN).LT.1.D-7) GO TO 3050                       03490000
      ASCA=0.D0                                                         03500000
      DO 7000 NA=1,NMAX                                                 03510000
      DO 7010 NB=1,NMAX                                                 03520000
      DSCA=ALP(NA)*DCONJG(ALP(NB))+BET(NA)*DCONJG(BET(NB))              03530000
      GSCA=BSCA(NA,NB)*DSCA                                             03540000
      ASCA=ASCA+GSCA                                                    03550000
 7010 CONTINUE                                                          03560000
 7000 CONTINUE                                                          03570000
      BEXT=0.D0                                                         03580000
      DO 7070 NA=1,NMAX                                                 03590000
      CEXT=ALP(NA)*SIG0(NA,IN)+BET(NA)*CHI0(NA,IN)                      03600000
      BEXT=BEXT+CEXT                                                    03610000
 7070 CONTINUE                                                          03620000
 7080 BEXT=-4.D0*BEXT                                                   03630000
      EXT(IN,IMODE)=EXT(IN,IMODE)+BEXT                                  03640000
      SCA(IN,IMODE)=SCA(IN,IMODE)+ASCA                                  03650000
      CONV(IN)=DABS(BEXT/EXT(IN,IMODE))                                 03660000
      IF(IMODE.EQ.2) GO TO 3050                                         03670000
      IF(IN.GT.1.AND.IN.LT.INCANG) GO TO 3050                           03680000
      WRITE(6,101) IN,DZETA(IN),BEXT,ASCA,NMAX                          03690000
  101 FORMAT(1H ,14X,1H(,I2,21H )   INCIDENT ANGLE =,F5.1,5H DEG. ,     03700000
     1 10X,6HAEXT =,1PD12.5,7X,6HASCA =,D12.5,5X,7H(NMAX =,I3,2H ) )    03710000
 3050 CONTINUE                                                          03720000
 5000 CONTINUE                                                          03730000
      CRIT1=EXT(IZETA,1)+EXT(IZETA,2)                                   03740000
      CRIT=DABS((CRIT1-CRIT0)/CRIT1)                                    03750000
      IF(M.GE.MMAX) GO TO 7                                             03760000
      IF(CRIT.LE.0.0001) GO TO 7                                        03770000
      M=M+1                                                             03780000
      GO TO 9050                                                        03790000
    7 MFIN=MM                                                           03800000
       REWIND NUNIT                                                     03810000
      DO 8000 IN=1,INCANG                                               03820000
      CRITM=1.D0                                                        03830000
CC            PRINT OUT OF THE CROSS SECTION AND EFFICIENCY FACTOR      03840000
      WRITE(6,600) DZETA(IN)                                            03850000
      XG=DSQRT(PG(IN)/PI)                                               03860000
      DO 8002 MODE=1,2                                                  03870000
      CEXT=PI*EXT(IN,MODE)                                              03880000
      CSCA=PI*SCA(IN,MODE)                                              03890000
      QEXT=CEXT/PG(IN)                                                  03900000
      QSCA=CSCA/PG(IN)                                                  03910000
      QABS=QEXT-QSCA                                                    03920000
      ALBEDO=QSCA/QEXT                                                  03930000
      IF(MODE.EQ.2) GO TO 8004                                          03940000
      WRITE(6,601)                                                      03950000
      GO TO 8006                                                        03960000
 8004 WRITE(6,602)                                                      03970000
 8006 WRITE(6,603) CEXT,CSCA,PG(IN),XG                                  03980000
      WRITE(6,604) QEXT,QSCA,QABS,ALBEDO                                03990000
 8002 CONTINUE                                                          04000000
      AVEXT=0.5*PI*(EXT(IN,1)+EXT(IN,2))                                04010000
      AVSCA=0.5*PI*(SCA(IN,1)+SCA(IN,2))                                04020000
      AQEXT=AVEXT/PG(IN)                                                04030000
      AQSCA=AVSCA/PG(IN)                                                04040000
      AQABS=AQEXT-AQSCA                                                 04050000
      AVALB=AQSCA/AQEXT                                                 04060000
      WRITE(6,609)                                                      04070000
      WRITE(6,603) AVEXT,AVSCA,PG(IN),XG                                04080000
      WRITE(6,604) AQEXT,AQSCA,AQABS,AVALB                              04090000
  600 FORMAT(1H0//10X,40('-'), ' INCIDENCE ANGLE =',F6.2,' DEG.',40('-')04100000
     + //12X,'*   CROSS SECTIONS AND EFFICIENCY FACTORS FOR EXTINCTION  04110000
     + AND SCATTERING' )                                                04120000
  601 FORMAT(1H /15X,'A)   TM-MODE POLARIZATION')                       04130000
  602 FORMAT(1H /15X,'B)   TE-MODE POLARIZATION')                       04140000
  603 FORMAT(1H ,15X,'CEXT =',1PD13.5,10X,'CSCA =',D13.5,10X,'AREA =',  04150000
     + D13.5,10X,'XG   =',D13.5 )                                       04160000
  604 FORMAT(1H ,15X,'QEXT =',1PD13.5,10X,'QSCA =',D13.5,10X,'QABS =',  04170000
     + D13.5,10X,'ALBD =',D13.5 )                                       04180000
  609 FORMAT(1H /15X,'C)   POLARIZATION AVERAGED CROSS SECTIONS' )      04190000
CC            AMPLITUDE FUNCTIONS AND TRANSFORMATION MATRIXES           04200000
C                                                                       04210031
        IF(ICALC.LT.0)  GO TO 8000                                      04220031
C                                                                       04230031
      INDEX1=IN                                                         04240031
      DO 2060 MM=1,MFIN                                                 04250000
      M=MM-1                                                            04260000
      NMAX=NUMB(MM)                                                     04270000
      DO 2060 MODE=1,2                                                  04280000
       READ(NUNIT'INDEX1) DATA                                          04290000
      DO 2070 NC=1,NMAX                                                 04300000
      ALP1(NC,MM,MODE)=DATA(NC)                                         04310000
      BET1(NC,MM,MODE)=DATA(NC+NMX)                                     04320000
 2070 CONTINUE                                                          04330000
      INDEX1=INDEX1+INCANG                                              04340000
 2060 CONTINUE                                                          04350000
      DO 8030 KP=1,NPHI1                                                04360000
      DO 8008 JT=1,NTHT                                                 04370000
        CALL ANGLE(ZETA(IN),0.D0,RLTHT(JT),RLPHI(KP),STHT,SPHI,SALPH)   04380000
      S2AL=DSIN(2.*SALPH)/2.                                            04390000
      CCAL=DCOS(SALPH)**2                                               04400000
      SSAL=DSIN(SALPH)**2                                               04410000
      A1=(0.D0,0.D0)                                                    04420000
      A2=(0.D0,0.D0)                                                    04430000
      A3=(0.D0,0.D0)                                                    04440000
      A4=(0.D0,0.D0)                                                    04450000
      ETA=DCOS(STHT)                                                    04460000
      ETA2=1.D0-ETA*ETA                                                 04470000
      SQRE=DSQRT(ETA2)                                                  04480000
      ABSETA=DABS(ETA)                                                  04490000
      DO 8070 MM=1,MFIN                                                 04500000
      M=MM-1                                                            04510000
      NMAX=NUMB(MM)                                                     04520000
      TE1=(0.D0,0.D0)                                                   04530000
      TE2=(0.D0,0.D0)                                                   04540000
      TM1=(0.D0,0.D0)                                                   04550000
      TM2=(0.D0,0.D0)                                                   04560000
      IF(ZETA(IN).LE.1.D-7.AND.M.NE.1) GO TO 8020                       04570000
      IF(M.GE.3.AND.CRITM.LT.1.D-6) GO TO 8020                          04580000
      IF(ABSETA.GE.OMG) GO TO 4000                                      04590000
      PL(1)=PFACT(MM)*(SQRE**M)                                         04600000
      PL(2)=(2*M+1)*ETA*PL(1)                                           04610000
      DO 4010 I=1,MAXJ1                                                 04620000
      MI=M+I                                                            04630000
      FI1=I+1                                                           04640000
      PL(I+2)=((2*MI+1)*ETA*PL(I+1)-(MI+M)*PL(I))/FI1                   04650000
 4010 DPL(I)=(MI*ETA*PL(I)-(MI-M)*PL(I+1))/ETA2                         04660000
 4000 CONTINUE                                                          04670000
      DO 8010 NC=1,NMAX                                                 04680000
      N=M+NC-1                                                          04690000
      NMMOD=MOD(N-M,2)                                                  04700000
      S=0.D0                                                            04710000
      DS=0.D0                                                           04720000
      SUM=0.D0                                                          04730000
      JEN=JWM(MM,NC)                                                    04740000
      DO 8040 IR=1,JEN,2                                                04750000
      I=IR                                                              04760000
      IF(NMMOD.EQ.0) I=IR-1                                             04770000
      IP=(I/2)+1                                                        04780000
      IF(ABSETA.GE.OMG) GO TO 8060                                      04790000
      DDPL=DEM(MM,NC,IP)*DPL(I+1)                                       04800000
      DS=DS+DDPL                                                        04810000
      IF(I.GT.NC.AND.DABS(DDPL).LT.1.D-30) GO TO 8050                   04820000
      IF(M.EQ.0) GO TO 8040                                             04830000
      S=S+DEM(MM,NC,IP)*PL(I+1)                                         04840000
      GO TO 8040                                                        04850000
 8060 SUM=SUM+.5D0*(I+1)*(I+2)*DEM(MM,NC,IP)                            04860000
 8040 CONTINUE                                                          04870000
 8050 IF(ABSETA.GE.OMG) GO TO 3                                         04880000
      SIG=M*S/SQRE                                                      04890000
      CHI=-SQRE*DS                                                      04900000
      GO TO 4                                                           04910000
    3 IF(M.NE.1) GO TO 8010                                             04920000
      IF(DABS(1.D0+ETA).LT.EPS) GO TO 6                                 04930000
      SIG=SUM                                                           04940000
      CHI=SIG                                                           04950000
      GO TO 4                                                           04960000
    6 SIG=SUM                                                           04970000
      IF(MOD(N-1,2).EQ.1) SIG=-SUM                                      04980000
      CHI=-SIG                                                          04990000
    4 TM1=TM1+ALP1(NC,MM,1)*CHI+BET1(NC,MM,1)*SIG                       05000000
      TM2=TM2+ALP1(NC,MM,1)*SIG+BET1(NC,MM,1)*CHI                       05010000
      TE1=TE1+ALP1(NC,MM,2)*SIG+BET1(NC,MM,2)*CHI                       05020000
      TE2=TE2+ALP1(NC,MM,2)*CHI+BET1(NC,MM,2)*SIG                       05030000
 8010 CONTINUE                                                          05040000
 8020 CONTINUE                                                          05050000
      FAI=M*SPHI                                                        05060000
      SINM=DSIN(FAI)                                                    05070000
      COSM=DCOS(FAI)                                                    05080000
      IF(DABS(SINM).LT.EPS) SINM=0.D0                                   05090000
      IF(DABS(COSM).LT.EPS) COSM=0.D0                                   05100000
      A1=A1+TE1*COSM                                                    05110000
      A2=A2+TM2*COSM                                                    05120000
      A3=A3-TE2*SINM                                                    05130000
      A4=A4+TM1*SINM                                                    05140000
      IF(M.EQ.0) GO TO 8070                                             05150000
      CRITM=CDABS(TM1+TM2)/CDABS(A4+A2)                                 05160000
 8070 CONTINUE                                                          05170000
CC        TRANSFORMATION MATRIX - F' FOR THE STOKES PARAM. (IL,IR)      05180000
      AR1= DREAL(A1)                                                    05190000
      AR2= DREAL(A2)                                                    05200000
      AR3= DREAL(A3)                                                    05210000
      AR4= DREAL(A4)                                                    05220000
      AI1= DIMAG(A1)                                                    05230000
      AI2= DIMAG(A2)                                                    05240000
      AI3= DIMAG(A3)                                                    05250000
      AI4= DIMAG(A4)                                                    05260000
      FD11=(AR2*AR2+AI2*AI2)                                            05270000
      FD12=(AR3*AR3+AI3*AI3)                                            05280000
      FD21=(AR4*AR4+AI4*AI4)                                            05290000
      FD22=(AR1*AR1+AI1*AI1)                                            05300000
      FD31=2.D0*(AR2*AR4+AI2*AI4)                                       05310000
      FD32=2.D0*(AR3*AR1+AI3*AI1)                                       05320000
C           INTENSITY MATRIX 'R' FOR (IL,IR) AT (THETA,PHI)             05330000
      R(JT,KP,1,1)=CCAL*FD11+SSAL*FD21+S2AL*FD31                        05340000
      R(JT,KP,1,2)=CCAL*FD12+SSAL*FD22+S2AL*FD32                        05350000
      R(JT,KP,2,1)=SSAL*FD11+CCAL*FD21-S2AL*FD31                        05360000
      R(JT,KP,2,2)=SSAL*FD12+CCAL*FD22-S2AL*FD32                        05370000
 8008 CONTINUE                                                          05380000
 8030 CONTINUE                                                          05390000
CC          AZIMUTHAL ANISOTROPY OF SCATTERED FIELDS                    05400000
      WRITE(6,605)                                                      05410000
      DO 4060 K=1,5                                                     05420000
      KK=NAZM(K)                                                        05430000
      AZMTH=RLPHI(KK)/STRD                                              05440000
      WRITE(6,606) K,AZMTH                                              05450000
      DO 4060 J=1,NTHT                                                  05460000
      TOTAL=(R(J,KK,1,1)+R(J,KK,1,2)+R(J,KK,2,1)+R(J,KK,2,2))/2.D0      05470000
      DGPOL=0.5*(R(J,KK,2,1)+R(J,KK,2,2)-R(J,KK,1,1)-R(J,KK,1,2))/TOTAL 05480000
      IF(R(J,KK,1,1).LT.1.D-20) GO TO 4040                              05490000
      DEPHV=R(J,KK,2,1)/R(J,KK,1,1)                                     05500000
      DEPVH=R(J,KK,1,2)/R(J,KK,2,2)                                     05510000
      GO TO 4050                                                        05520000
 4040 DEPHV=1./EPS                                                      05530000
      DEPVH=1./EPS                                                      05540000
 4050 WRITE(6,607) SANG(J),((R(J,KK,I,L),L=1,2),I=1,2),TOTAL,DGPOL,     05550000
     + DEPHV,DEPVH                                                      05560000
  605 FORMAT(1H0/10X,'*   AZIMUTHAL ANISOTROPY OF THE SCATTERED INTENSIT05570000
     +Y AND DEGREE OF POLARIZATION FOR NATURAL LIGHT' )                 05580000
  606 FORMAT(1H /15X,I1,')   AZIMUTH ANGLE =',F7.2,' DEG.'/4X,'SANG',   05590000
     + 8X,'I-LL',11X,'I-RL',11X,'I-LR',11X,'I-RR',11X,'TOTAL',10X,      05600000
     + 'DGPOL',8X,'I-LR/I-LL',6X,'I-RL/I-RR' )                          05610000
  607 FORMAT(1H ,F7.2,1P8D15.4)                                         05620000
 4060 CONTINUE                                                          05630000
CC        AZIMUTHALLY AVERAGED INTENSITY MATRIX 'P'                     05640000
      DO 4070 JT=1,NTHT                                                 05650000
      DO 4080 J=1,2                                                     05660000
      DO 4080 I=1,2                                                     05670000
      P(JT,I,J)=0.D0                                                    05680000
      DO 4090 KP=1,NPHI1                                                05690000
 4090 P(JT,I,J)=P(JT,I,J)+WPHI(KP)*R(JT,KP,I,J)                         05700000
 4080 CONTINUE                                                          05710000
 4070 CONTINUE                                                          05720000
CC        NORMALIZATION OF THE INTENSITY MATRIX                         05730000
      WRITE(6,801)                                                      05740000
      WRITE(6,802)                                                      05750000
  801 FORMAT(1H0/10X,   '*   AZIMUTHALLY AVERAGED, NORMALIZED INTENSITY 05760000
     +MATRIX - P ' )                                                    05770000
  802 FORMAT(1H0,4X,'SANG',6X,'P(1,1)',9X,'P(1,2)',9X,'P(2,1)',9X,      05780000
     1 'P(2,2)',10X,'TOTAL',10X,'DGPOL',9X,'P21/P11',8X,'P12/P22' )     05790000
      PNORM=0.D0                                                        05800000
      COSBAR=0.D0                                                       05810000
      CSCAT=PI*(SCA(IN,1)+SCA(IN,2))/2.                                 05820000
      CNORM=4.*PI/CSCAT                                                 05830000
      DO 2020 J=1,NTHT                                                  05840000
      WW=DSIN(RLTHT(J))*WTHT(J)                                         05850000
      DO 2030 M=1,2                                                     05860000
      DO 2030 N=1,2                                                     05870000
 2030 P(J,M,N)=CNORM*P(J,M,N)                                           05880000
      TOTAL=(P(J,1,1)+P(J,1,2)+P(J,2,1)+P(J,2,2))/2.                    05890000
      DGPOL=0.5*(P(J,2,1)+P(J,2,2)-P(J,1,2)-P(J,1,1))/TOTAL             05900000
      IF(P(J,1,1).LT.1.D-20) GO TO 2040                                 05910000
      DEPHV=P(J,2,1)/P(J,1,1)                                           05920000
      DEPVH=P(J,1,2)/P(J,2,2)                                           05930000
      GO TO 2050                                                        05940000
 2040 DEPHV=1./EPS                                                      05950000
      DEPVH=1./EPS                                                      05960000
 2050 PNORM=PNORM+WW*TOTAL                                              05970000
      COSBAR=COSBAR+WW*DCOS(RLTHT(J))*TOTAL                             05980000
      WRITE(6,607) SANG(J),((P(J,K,L),L=1,2),K=1,2),TOTAL,DGPOL,        05990000
     + DEPHV,DEPVH                                                      06000000
 2020 CONTINUE                                                          06010000
      COSBAR=COSBAR/PNORM                                               06020000
      PNORM=2.D0*PI*PNORM/(4.D0*PI)                                     06030000
      WRITE(6,803) COSBAR,PNORM                                         06040000
  803 FORMAT(1H0/10X,'*     ASYMMETRY PARAMETER  -  <COS> =',1PD13.5/   06050000
     1 16X,'NORMALIZATION =',D13.5 )                                    06060000
 8000 CONTINUE                                                          06070000
C      CALL CLOCK(CPTIME)                                               06080000
C      WRITE(6,6655) CPTIME                                             06090000
C6655 FORMAT(1H0/100X,'CPU TIME =', F10.3,'SEC' )                       06100000
      GO TO 400                                                         06110000
  300 STOP                                                              06120000
      END                                                               06130000

CCCCC  -----  Example of Input Data Cards  -------------------
 DINC= 45.0 DSANG=  1.0 DPHI= 10.0
 NSHAPE= 1 ABYB=  1.5 A=  2.0 RFRRL= 1.450 RFRIM= 0.0
 NSHAPE= 1 ABYB=  2.0 A= 10.0 RFRRL= 1.330 RFRIM= 1.0E-3
 NSHAPE= 2 ABYB=  1.5 A=  5.0 RFRRL= 1.450 RFRIM= 0.0
 NSHAPE= 2 ABYB=  2.0 A= 10.0 RFRRL= 1.330 RFRIM= 1.0E-3
CCCCC  ------  End of Example Input Data  --------------------  